knitr::opts_chunk$set(echo = TRUE)

library(tidyverse) # great collection of packages for data carpentry, modelling, and visualization
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.3     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(readr)
library(haven) # package for loading Stata's .dta files. 
library(sjlabelled) # good for renaming, changing classes, etc. in the piped dplyr mode
## 
## Attaching package: 'sjlabelled'
## The following objects are masked from 'package:haven':
## 
##     as_factor, read_sas, read_spss, read_stata, write_sas, zap_labels
## The following object is masked from 'package:forcats':
## 
##     as_factor
## The following object is masked from 'package:dplyr':
## 
##     as_label
## The following object is masked from 'package:ggplot2':
## 
##     as_label
library(zscorer)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(anthro)
library(childsds)
library(sjlabelled)
library(sjPlot)
## Registered S3 methods overwritten by 'parameters':
##   method                           from      
##   as.double.parameters_kurtosis    datawizard
##   as.double.parameters_skewness    datawizard
##   as.double.parameters_smoothness  datawizard
##   as.numeric.parameters_kurtosis   datawizard
##   as.numeric.parameters_skewness   datawizard
##   as.numeric.parameters_smoothness datawizard
##   print.parameters_distribution    datawizard
##   print.parameters_kurtosis        datawizard
##   print.parameters_skewness        datawizard
##   summary.parameters_kurtosis      datawizard
##   summary.parameters_skewness      datawizard
library(VIF)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:VIF':
## 
##     vif
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
growth_clocks_data <- read_csv(here::here ("Data/growth_clocks_data.csv"))
## New names:
## * `` -> ...1
## Rows: 3023 Columns: 158
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr    (9): basebrgy_basewman, reprostat, was_preg_no_na, trimester, SampleI...
## dbl  (144): ...1, uncchdid, sexchild, sex, outcome, age_days_birthweigh, wei...
## date   (5): birthdate, intwdate_birth, date_inf_meas, intwdate91, intw_date02
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
growth_clocks_data
## # A tibble: 3,023 × 158
##     ...1 uncchdid basebrgy_basewman sexchild   sex outcome birthdate 
##    <dbl>    <dbl> <chr>                <dbl> <dbl>   <dbl> <date>    
##  1     1    20001 1_12                     1    NA       1 1983-06-03
##  2     2    20002 1_14                     2    NA       1 1983-05-23
##  3     3    20003 1_21                     2    NA       1 1983-05-13
##  4     4    20004 1_23                     1     1       1 1983-05-19
##  5     5    20006 1_26                     2    NA       1 1983-06-20
##  6     6    20007 1_27                     2     2       1 1983-05-16
##  7     7    20008 1_30                     1     1       1 1983-05-03
##  8     8    20009 1_32                     2    NA       1 1983-05-20
##  9     9    20010 1_33                     2     2       1 1983-05-07
## 10    10    20011 1_34                     2    NA       1 1983-05-31
## # … with 3,013 more rows, and 151 more variables: intwdate_birth <date>,
## #   age_days_birthweigh <dbl>, weightak_birth_kg <dbl>, heightcm_birth <dbl>,
## #   date_inf_meas <date>, age_days_infweigh <dbl>, hght_12 <dbl>,
## #   wght_12 <dbl>, intwdate91 <date>, age91_days <dbl>, age_mo_91 <dbl>,
## #   age91_years <dbl>, hightndx_91 <dbl>, weghtndx_91 <dbl>,
## #   intw_date02 <date>, age_days_02 <dbl>, age_years_02 <dbl>,
## #   age_cutoff_02 <dbl>, age_cutoff_02_days <dbl>, height_02 <dbl>, …

#Birth to 2 years old

grim_height_b_2_f<-lm(AgeAccelGrim ~ hfa_diff_birth_inf12 + 
                        was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_height_b_2_m<-lm(AgeAccelGrim ~ hfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_height_b_2_f <-update(grim_height_b_2_f, AgeAccelPheno ~ .)

pheno_height_b_2_m <-update(grim_height_b_2_m, AgeAccelPheno ~ .-was_preg_no_na)

han_height_b_2_f <-update(grim_height_b_2_f, EEAA ~ .)


han_height_b_2_m <-update(grim_height_b_2_m, EEAA ~ . -was_preg_no_na)


horv_height_b_2_f <-update(grim_height_b_2_f, IEAA ~ .)


horv_height_b_2_m <-update(grim_height_b_2_m, IEAA ~ . -was_preg_no_na)

sjPlot::tab_model(grim_height_b_2_f, grim_height_b_2_m, 
                  pheno_height_b_2_f, pheno_height_b_2_m, 
                  han_height_b_2_f, han_height_b_2_m, 
                  horv_height_b_2_f, horv_height_b_2_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.93 -1.39 – -0.48 <0.001 1.80 0.76 – 2.84 0.001 -0.34 -1.29 – 0.61 0.478 -2.19 -4.02 – -0.37 0.019 -0.39 -1.39 – 0.61 0.442 0.33 -1.49 – 2.16 0.719 0.09 -0.50 – 0.68 0.757 1.15 -0.16 – 2.47 0.085
hfa_diff_birth_inf12 -0.04 -0.23 – 0.16 0.717 0.16 -0.30 – 0.63 0.483 -0.04 -0.43 – 0.36 0.859 -0.51 -1.33 – 0.30 0.214 0.01 -0.41 – 0.43 0.962 -0.35 -1.17 – 0.47 0.397 0.09 -0.16 – 0.34 0.473 0.54 -0.04 – 1.13 0.070
was_preg_no_na [Yes] 2.79 2.20 – 3.39 <0.001 3.47 2.23 – 4.71 <0.001 0.96 -0.34 – 2.26 0.146 -0.00 -0.77 – 0.77 0.995
Observations 372 100 372 100 372 100 372 100
R2 / R2 adjusted 0.189 / 0.185 0.005 / -0.005 0.076 / 0.071 0.016 / 0.006 0.006 / 0.000 0.007 / -0.003 0.001 / -0.004 0.033 / 0.023 
growth_clocks_data %>% 
  select(uncchdid, hfa_diff_birth_inf12, AgeAccelGrim, AgeAccelPheno, EEAA, IEAA) %>% 
  gather(key = clock_type, value = AgeAccel, -c(1,2)) %>% 
  na.omit() %>% 
  ggplot(., aes(x = hfa_diff_birth_inf12, y = AgeAccel, col = clock_type))+
  geom_point()+
  scale_color_brewer(type = "qual", palette = 6)+
  facet_wrap(~clock_type)+
  theme(legend.position = "none")

#birth-2 yrs old visualization 
par(mfrow=c(2,2))
plot(grim_height_b_2_f)

par(mfrow=c(2,2))
plot(grim_height_b_2_m)

par(mfrow=c(2,2))
plot(pheno_height_b_2_f)

par(mfrow=c(2,2))
plot(pheno_height_b_2_m)

par(mfrow=c(2,2))
plot(han_height_b_2_f)

par(mfrow=c(2,2))
plot(han_height_b_2_m)

#hfaz 2years old model
grim_height_83_91_f <-lm(AgeAccelGrim ~ hfa_diff_inf12_91 + 
                        hfa_diff_birth_inf12 +
                        was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_height_83_91_m<-lm(AgeAccelGrim ~ hfa_diff_inf12_91 + 
                        hfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_height_83_91_f <-update(grim_height_83_91_f, AgeAccelPheno ~ .)

pheno_height_83_91_m <-update(grim_height_83_91_m, AgeAccelPheno ~ .-was_preg_no_na)

han_height_83_91_f <-update(grim_height_83_91_f, EEAA ~ .)

han_height_83_91_m <-update(grim_height_83_91_m, EEAA ~ .-was_preg_no_na)

horv_height_83_91_f <-update(grim_height_83_91_f, IEAA ~ .)

horv_height_83_91_m <-update(grim_height_83_91_m, IEAA ~ .-was_preg_no_na)


sjPlot::tab_model(grim_height_83_91_f, grim_height_83_91_m, 
                  pheno_height_83_91_f, pheno_height_83_91_m, 
                  han_height_83_91_f, han_height_83_91_m, 
                  horv_height_83_91_f, horv_height_83_91_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.89 -1.35 – -0.43 <0.001 1.96 0.93 – 2.98 <0.001 -0.24 -1.20 – 0.71 0.617 -1.86 -3.64 – -0.09 0.040 -0.21 -1.21 – 0.78 0.672 0.71 -1.05 – 2.46 0.426 0.07 -0.53 – 0.66 0.820 1.28 -0.04 – 2.60 0.056
hfa_diff_inf12_91 0.16 -0.12 – 0.43 0.268 0.92 0.16 – 1.68 0.018 0.18 -0.40 – 0.76 0.542 1.93 0.61 – 3.24 0.004 0.52 -0.09 – 1.12 0.093 2.19 0.89 – 3.48 0.001 -0.31 -0.67 – 0.05 0.094 0.75 -0.23 – 1.73 0.131
hfa_diff_birth_inf12 0.01 -0.19 – 0.22 0.898 0.41 -0.08 – 0.91 0.103 0.03 -0.40 – 0.47 0.876 0.01 -0.86 – 0.87 0.989 0.18 -0.27 – 0.63 0.433 0.24 -0.61 – 1.09 0.576 0.01 -0.26 – 0.28 0.953 0.75 0.11 – 1.39 0.023
was_preg_no_na [Yes] 2.76 2.17 – 3.36 <0.001 3.42 2.18 – 4.66 <0.001 0.85 -0.44 – 2.14 0.196 0.02 -0.75 – 0.79 0.957
Observations 370 100 370 100 370 100 370 100
R2 / R2 adjusted 0.191 / 0.185 0.061 / 0.042 0.076 / 0.069 0.095 / 0.076 0.013 / 0.005 0.110 / 0.092 0.009 / 0.001 0.056 / 0.036 
#83-91 visualization 
par(mfrow=c(2,2))
plot(grim_height_83_91_f)

par(mfrow=c(2,2))
plot(grim_height_83_91_m)

par(mfrow=c(2,2))
plot(pheno_height_83_91_f)

par(mfrow=c(2,2))
plot(pheno_height_83_91_m)

par(mfrow=c(2,2))
plot(han_height_83_91_f)

par(mfrow=c(2,2))
plot(han_height_83_91_m)

8 to 19 years old

#hfaz_91 minimal models

grim_height_91_02_f<-lm(AgeAccelGrim ~ hfa_diff_91_02 +
                     hfa_diff_birth_inf12 +
                     was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_height_91_02_m<-lm(AgeAccelGrim ~ hfa_diff_91_02 +
                     hfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_height_91_02_f <-update(grim_height_91_02_f, AgeAccelPheno ~ .)

pheno_height_91_02_m <-update(grim_height_91_02_m, AgeAccelPheno ~ .-was_preg_no_na)

han_height_91_02_f <-update(grim_height_91_02_f, EEAA ~ .)

han_height_91_02_m <-update(grim_height_91_02_m, EEAA ~ .-was_preg_no_na)

horv_height_91_02_f <-update(grim_height_91_02_f, IEAA ~ .)

horv_height_91_02_m <-update(grim_height_91_02_m, IEAA ~ .-was_preg_no_na)


sjPlot::tab_model(grim_height_91_02_f, grim_height_91_02_m, 
                  pheno_height_91_02_f, pheno_height_91_02_m, 
                  han_height_91_02_f, han_height_91_02_m, 
                  horv_height_91_02_f, horv_height_91_02_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.94 -1.39 – -0.48 <0.001 1.90 0.86 – 2.93 <0.001 -0.30 -1.26 – 0.66 0.538 -2.05 -3.87 – -0.22 0.029 -0.35 -1.35 – 0.65 0.493 0.33 -1.51 – 2.18 0.721 0.11 -0.47 – 0.70 0.703 1.21 -0.11 – 2.54 0.073
hfa_diff_91_02 -0.33 -0.68 – 0.03 0.069 -0.90 -1.97 – 0.16 0.096 -0.38 -1.13 – 0.36 0.310 -1.35 -3.22 – 0.53 0.158 -0.51 -1.29 – 0.26 0.194 -0.02 -1.91 – 1.88 0.987 0.28 -0.18 – 0.74 0.229 -0.54 -1.90 – 0.83 0.435
hfa_diff_birth_inf12 -0.07 -0.27 – 0.13 0.480 0.11 -0.36 – 0.57 0.644 -0.05 -0.46 – 0.36 0.820 -0.60 -1.42 – 0.22 0.151 -0.04 -0.47 – 0.38 0.840 -0.35 -1.18 – 0.48 0.403 0.15 -0.10 – 0.41 0.238 0.51 -0.08 – 1.11 0.092
was_preg_no_na [Yes] 2.72 2.12 – 3.32 <0.001 3.45 2.20 – 4.70 <0.001 0.87 -0.44 – 2.17 0.192 0.07 -0.70 – 0.84 0.859
Observations 364 100 364 100 364 100 364 100
R2 / R2 adjusted 0.193 / 0.187 0.033 / 0.013 0.080 / 0.073 0.036 / 0.016 0.010 / 0.002 0.007 / -0.013 0.006 / -0.002 0.039 / 0.019 
#91-02 visualization
par(mfrow=c(2,2))
plot(grim_height_91_02_f)

par(mfrow=c(2,2))
plot(grim_height_91_02_m)

par(mfrow=c(2,2))
plot(pheno_height_91_02_f)

par(mfrow=c(2,2))
plot(pheno_height_91_02_m)

par(mfrow=c(2,2))
plot(han_height_91_02_f)

par(mfrow=c(2,2))
plot(han_height_91_02_m)

#modeling wfaz (no interactions)

Birth to 2 years old

#wfa birth minimal models
grim_weight_b_2_f<-lm(AgeAccelGrim ~ wfa_diff_birth_inf12 + was_preg_no_na, subset(growth_clocks_data, sex == "2")) 

grim_weight_b_2_m<-lm(AgeAccelGrim ~ wfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_weight_b_2_f <-update(grim_weight_b_2_f, AgeAccelPheno ~ .)

pheno_weight_b_2_m <-update(grim_weight_b_2_m, AgeAccelPheno ~ .)

han_weight_b_2_f <-update(grim_weight_b_2_f, EEAA ~ .)

han_weight_b_2_m <-update(grim_weight_b_2_m, EEAA ~ .)

horv_weight_b_2_f <-update(grim_weight_b_2_f, IEAA ~ .)

horv_weight_b_2_m <-update(grim_weight_b_2_m, IEAA ~ .)


sjPlot::tab_model(grim_weight_b_2_f, grim_weight_b_2_m, 
                  pheno_weight_b_2_f, pheno_weight_b_2_m, 
                  han_weight_b_2_f, han_weight_b_2_m, 
                  horv_weight_b_2_f, horv_weight_b_2_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.76 -1.46 – -0.07 0.031 1.10 -0.84 – 3.04 0.262 0.29 -1.17 – 1.74 0.698 0.32 -3.09 – 3.73 0.851 0.76 -0.76 – 2.28 0.327 1.90 -1.51 – 5.31 0.271 0.30 -0.61 – 1.20 0.521 -0.37 -2.86 – 2.12 0.771
wfa_diff_birth_inf12 -0.04 -0.30 – 0.22 0.770 0.14 -0.53 – 0.81 0.683 -0.22 -0.76 – 0.32 0.420 -0.56 -1.73 – 0.62 0.352 -0.46 -1.03 – 0.10 0.106 -0.33 -1.50 – 0.85 0.584 -0.15 -0.49 – 0.18 0.373 0.18 -0.68 – 1.04 0.684
was_preg_no_na [Yes] 2.80 2.21 – 3.40 <0.001 3.52 2.28 – 4.76 <0.001 1.05 -0.25 – 2.35 0.112 0.03 -0.75 – 0.80 0.948
Observations 372 100 372 100 372 100 372 100
R2 / R2 adjusted 0.189 / 0.185 0.002 / -0.008 0.078 / 0.073 0.009 / -0.001 0.013 / 0.007 0.003 / -0.007 0.002 / -0.003 0.002 / -0.008 
#weight birth-2 yrs old 
par(mfrow=c(2,2))
plot(grim_weight_b_2_f)

par(mfrow=c(2,2))
plot(grim_weight_b_2_m)

par(mfrow=c(2,2))
plot(pheno_weight_b_2_f)

par(mfrow=c(2,2))
plot(pheno_weight_b_2_m)

par(mfrow=c(2,2))
plot(han_weight_b_2_f)

par(mfrow=c(2,2))
plot(han_weight_b_2_m)

2 to 8 years old

#wfaz inf12 minimal models
grim_weight_83_91_f<-lm(AgeAccelGrim ~ wfa_diff_inf12_91 + 
                          wfa_diff_birth_inf12+
                          was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_weight_83_91_m<-lm(AgeAccelGrim ~ wfa_diff_inf12_91+
                          wfa_diff_birth_inf12,
                          subset(growth_clocks_data, sex == "1"))
                        
pheno_weight_83_91_f <-update(grim_weight_83_91_f, AgeAccelPheno ~ .)

pheno_weight_83_91_m <-update(grim_weight_83_91_m, AgeAccelPheno ~ .)

han_weight_83_91_f <-update(grim_weight_83_91_f, EEAA ~ .)

han_weight_83_91_m <-update(grim_weight_83_91_m, EEAA ~ .)

horv_weight_83_91_f <-update(grim_weight_83_91_f, IEAA ~ .)

horv_weight_83_91_m <-update(grim_weight_83_91_m, IEAA ~ .)


sjPlot::tab_model(grim_weight_83_91_f, grim_weight_83_91_m, 
                  pheno_weight_83_91_f, pheno_weight_83_91_m, 
                  han_weight_83_91_f, han_weight_83_91_m, 
                  horv_weight_83_91_f, horv_weight_83_91_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.92 -1.68 – -0.16 0.018 1.02 -0.91 – 2.96 0.296 0.13 -1.46 – 1.71 0.875 0.37 -3.05 – 3.80 0.829 0.53 -1.12 – 2.18 0.532 1.83 -1.59 – 5.26 0.290 0.46 -0.53 – 1.45 0.361 -0.41 -2.91 – 2.10 0.748
wfa_diff_inf12_91 0.20 -0.12 – 0.52 0.219 0.58 -0.26 – 1.42 0.175 0.28 -0.39 – 0.95 0.412 -0.38 -1.88 – 1.11 0.612 0.41 -0.29 – 1.11 0.247 0.49 -1.00 – 1.98 0.520 -0.13 -0.55 – 0.29 0.537 0.29 -0.80 – 1.38 0.594
wfa_diff_birth_inf12 0.04 -0.26 – 0.33 0.803 0.24 -0.44 – 0.93 0.480 -0.13 -0.74 – 0.48 0.671 -0.63 -1.84 – 0.59 0.309 -0.33 -0.97 – 0.30 0.303 -0.24 -1.45 – 0.98 0.699 -0.22 -0.60 – 0.16 0.258 0.23 -0.66 – 1.12 0.607
was_preg_no_na [Yes] 2.78 2.18 – 3.37 <0.001 3.47 2.23 – 4.71 <0.001 0.98 -0.31 – 2.27 0.136 0.02 -0.75 – 0.80 0.954
Observations 370 100 370 100 370 100 370 100
R2 / R2 adjusted 0.192 / 0.185 0.021 / 0.000 0.079 / 0.072 0.011 / -0.009 0.017 / 0.009 0.007 / -0.013 0.004 / -0.005 0.005 / -0.016 
#weight 83-91 yrs old 
par(mfrow=c(2,2))
plot(grim_weight_83_91_f)

par(mfrow=c(2,2))
plot(grim_weight_83_91_m)

par(mfrow=c(2,2))
plot(pheno_weight_83_91_f)

par(mfrow=c(2,2))
plot(pheno_weight_83_91_m)

par(mfrow=c(2,2))
plot(han_weight_83_91_f)

par(mfrow=c(2,2))
plot(han_weight_83_91_m)

modeling hfaz (interactions)

2 to 8 years old

#hfaz 2years old model
grim_height_83_91_intxn_f<-lm(AgeAccelGrim ~ hfa_diff_inf12_91 * 
                        hfa_diff_birth_inf12 +
                        was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_height_83_91_intxn_m<-lm(AgeAccelGrim ~ hfa_diff_inf12_91 * 
                        hfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_height_83_91_intxn_f <-update(grim_height_83_91_intxn_f, AgeAccelPheno ~ .)

pheno_height_83_91_intxn_m <-update(grim_height_83_91_intxn_m, AgeAccelPheno ~ .-was_preg_no_na)

han_height_83_91_intxn_f <-update(grim_height_83_91_intxn_f, EEAA ~ .)

han_height_83_91_intxn_m <-update(grim_height_83_91_intxn_m, EEAA ~ .-was_preg_no_na)

horv_height_83_91_intxn_f <-update(grim_height_83_91_intxn_f, IEAA ~ .)

horv_height_83_91_intxn_m <-update(grim_height_83_91_intxn_m, IEAA ~ .-was_preg_no_na)


sjPlot::tab_model(grim_height_83_91_intxn_f, grim_height_83_91_intxn_m, 
                  pheno_height_83_91_intxn_f, pheno_height_83_91_intxn_m, 
                  han_height_83_91_intxn_f, han_height_83_91_intxn_m, 
                  horv_height_83_91_intxn_f, horv_height_83_91_intxn_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.92 -1.41 – -0.43 <0.001 1.75 0.71 – 2.80 0.001 -0.18 -1.20 – 0.84 0.730 -1.81 -3.65 – 0.02 0.052 -0.11 -1.17 – 0.96 0.844 0.60 -1.21 – 2.41 0.512 -0.04 -0.67 – 0.60 0.904 1.03 -0.31 – 2.38 0.131
hfa_diff_inf12_91 0.25 -0.27 – 0.76 0.344 1.72 0.50 – 2.93 0.006 0.02 -1.06 – 1.09 0.975 1.73 -0.40 – 3.86 0.110 0.24 -0.88 – 1.36 0.671 2.60 0.50 – 4.71 0.016 -0.03 -0.70 – 0.63 0.919 1.73 0.17 – 3.30 0.031
hfa_diff_birth_inf12 -0.01 -0.25 – 0.23 0.926 0.22 -0.33 – 0.76 0.428 0.08 -0.42 – 0.58 0.756 0.05 -0.90 – 1.01 0.912 0.25 -0.26 – 0.77 0.335 0.14 -0.81 – 1.08 0.771 -0.07 -0.38 – 0.24 0.673 0.51 -0.20 – 1.21 0.156
was_preg_no_na [Yes] 2.76 2.17 – 3.35 <0.001 3.42 2.18 – 4.66 <0.001 0.86 -0.43 – 2.15 0.192 0.01 -0.76 – 0.78 0.974
hfa_diff_inf12_91 *
hfa_diff_birth_inf12 		0.04 -0.16 – 0.24 0.680 0.43 -0.08 – 0.94 0.100 -0.07 -0.49 – 0.34 0.723 -0.10 -1.00 – 0.79 0.819 -0.13 -0.55 – 0.30 0.566 0.22 -0.66 – 1.11 0.618 0.12 -0.13 – 0.38 0.339 0.52 -0.13 – 1.18 0.117
Observations 370 100 370 100 370 100 370 100
R2 / R2 adjusted 0.192 / 0.183 0.087 / 0.059 0.077 / 0.067 0.095 / 0.067 0.014 / 0.003 0.112 / 0.085 0.012 / 0.001 0.080 / 0.051 
#height 83-91 interaction models 
par(mfrow=c(2,2))
plot(grim_height_83_91_intxn_f)

par(mfrow=c(2,2))
plot(grim_height_83_91_intxn_m)

par(mfrow=c(2,2))
plot(pheno_height_83_91_intxn_f)

par(mfrow=c(2,2))
plot(pheno_height_83_91_intxn_m)

par(mfrow=c(2,2))
plot(han_height_83_91_intxn_f)

par(mfrow=c(2,2))
plot(han_height_83_91_intxn_m)

8 to 19 years old

#hfaz_91 minimal models

grim_height_91_02_intxn_f<-lm(AgeAccelGrim ~ hfa_diff_91_02 *
                     hfa_diff_birth_inf12 +
                     was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_height_91_02_intxn_m<-lm(AgeAccelGrim ~ hfa_diff_91_02 *
                     hfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_height_91_02_intxn_f <-update(grim_height_91_02_intxn_f, AgeAccelPheno ~ .)

pheno_height_91_02_intxn_m <-update(grim_height_91_02_intxn_m, AgeAccelPheno ~ .)

han_height_91_02_intxn_f <-update(grim_height_91_02_intxn_f, EEAA ~ .)

han_height_91_02_intxn_m <-update(grim_height_91_02_intxn_m, EEAA ~ .)

horv_height_91_02_intxn_f <-update(grim_height_91_02_intxn_f, IEAA ~ .)

horv_height_91_02_intxn_m <-update(grim_height_91_02_intxn_m, IEAA ~ .)

sjPlot::tab_model(grim_height_91_02_intxn_f, grim_height_91_02_intxn_m, 
                  pheno_height_91_02_intxn_f, pheno_height_91_02_intxn_m, 
                  han_height_91_02_intxn_f, han_height_91_02_intxn_m, 
                  horv_height_91_02_intxn_f, horv_height_91_02_intxn_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.89 -1.35 – -0.42 <0.001 2.08 1.00 – 3.16 <0.001 -0.22 -1.20 – 0.76 0.659 -1.57 -3.46 – 0.31 0.101 -0.39 -1.41 – 0.63 0.455 0.90 -1.00 – 2.80 0.349 0.19 -0.41 – 0.79 0.538 1.67 0.31 – 3.02 0.017
hfa_diff_91_02 -0.64 -1.34 – 0.06 0.071 -1.99 -4.12 – 0.15 0.068 -0.89 -2.35 – 0.58 0.233 -4.16 -7.90 – -0.43 0.029 -0.26 -1.79 – 1.27 0.737 -3.39 -7.14 – 0.36 0.076 -0.19 -1.09 – 0.71 0.676 -3.24 -5.92 – -0.55 0.019
hfa_diff_birth_inf12 -0.04 -0.24 – 0.17 0.740 0.24 -0.28 – 0.75 0.359 0.01 -0.43 – 0.44 0.966 -0.26 -1.16 – 0.64 0.571 -0.07 -0.53 – 0.38 0.754 0.06 -0.85 – 0.96 0.903 0.21 -0.06 – 0.47 0.132 0.84 0.19 – 1.48 0.012
was_preg_no_na [Yes] 2.70 2.11 – 3.30 <0.001 3.43 2.17 – 4.68 <0.001 0.88 -0.43 – 2.19 0.188 0.05 -0.72 – 0.82 0.901
hfa_diff_91_02 *
hfa_diff_birth_inf12 		-0.15 -0.44 – 0.14 0.304 -0.60 -1.63 – 0.43 0.248 -0.24 -0.85 – 0.36 0.432 -1.57 -3.37 – 0.23 0.087 0.12 -0.51 – 0.75 0.709 -1.87 -3.68 – -0.07 0.042 -0.23 -0.60 – 0.15 0.233 -1.50 -2.79 – -0.21 0.024
Observations 364 100 364 100 364 100 364 100
R2 / R2 adjusted 0.196 / 0.187 0.047 / 0.017 0.082 / 0.072 0.065 / 0.036 0.011 / -0.001 0.049 / 0.020 0.010 / -0.001 0.090 / 0.061 
#height 91-02 interaction models 
par(mfrow=c(2,2))
plot(grim_height_91_02_intxn_f)

par(mfrow=c(2,2))
plot(grim_height_91_02_intxn_m)

par(mfrow=c(2,2))
plot(pheno_height_91_02_intxn_f)

par(mfrow=c(2,2))
plot(pheno_height_91_02_intxn_m)

par(mfrow=c(2,2))
plot(han_height_91_02_intxn_f)

par(mfrow=c(2,2))
plot(han_height_91_02_intxn_m)

#modeling wfaz (interactions)

2 to 8 years old

#wfaz inf12 minimal models
grim_weight_83_91_intxn_f<-lm(AgeAccelGrim ~ wfa_diff_inf12_91 * wfa_diff_birth_inf12 + was_preg_no_na, subset(growth_clocks_data, sex == "2"))

grim_weight_83_91_intxn_m<-lm(AgeAccelGrim ~ wfa_diff_inf12_91 * wfa_diff_birth_inf12, subset(growth_clocks_data, sex == "1"))

pheno_weight_83_91_intxn_f <-update(grim_weight_83_91_intxn_f, AgeAccelPheno ~ .)

pheno_weight_83_91_intxn_m <-update(grim_weight_83_91_intxn_m, AgeAccelPheno ~ .)

han_weight_83_91_intxn_f <-update(grim_weight_83_91_intxn_f, EEAA ~ .)

han_weight_83_91_intxn_m <-update(grim_weight_83_91_intxn_m, EEAA ~ .)

horv_weight_83_91_intxn_f <-update(grim_weight_83_91_intxn_f, IEAA ~ .)

horv_weight_83_91_intxn_m <-update(grim_weight_83_91_intxn_m, IEAA ~ .)


sjPlot::tab_model(grim_weight_83_91_intxn_f, grim_weight_83_91_intxn_m, 
                  pheno_weight_83_91_intxn_f, pheno_weight_83_91_intxn_m, 
                  han_weight_83_91_intxn_f, han_weight_83_91_intxn_m, 
                  horv_weight_83_91_intxn_f, horv_weight_83_91_intxn_m)
  AgeAccelGrim AgeAccelGrim AgeAccelPheno AgeAccelPheno EEAA EEAA IEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -1.07 -1.85 – -0.28 0.008 0.95 -1.12 – 3.03 0.364 -0.26 -1.89 – 1.38 0.757 1.28 -2.36 – 4.93 0.486 0.21 -1.49 – 1.91 0.807 2.63 -1.01 – 6.28 0.155 0.50 -0.52 – 1.53 0.335 -0.60 -3.28 – 2.09 0.660
wfa_diff_inf12_91 0.66 -0.04 – 1.35 0.064 0.31 -2.59 – 3.21 0.834 1.47 0.03 – 2.92 0.046 3.09 -2.01 – 8.18 0.232 1.39 -0.12 – 2.89 0.072 3.54 -1.56 – 8.63 0.171 -0.26 -1.17 – 0.64 0.568 -0.44 -4.19 – 3.32 0.818
wfa_diff_birth_inf12 0.07 -0.23 – 0.36 0.650 0.27 -0.47 – 1.02 0.469 -0.05 -0.66 – 0.56 0.872 -0.98 -2.29 – 0.33 0.140 -0.27 -0.91 – 0.37 0.414 -0.55 -1.86 – 0.76 0.407 -0.23 -0.61 – 0.16 0.245 0.31 -0.66 – 1.27 0.531
was_preg_no_na [Yes] 2.79 2.19 – 3.38 <0.001 3.49 2.25 – 4.72 <0.001 1.00 -0.29 – 2.28 0.129 0.02 -0.76 – 0.80 0.958
wfa_diff_inf12_91 *
wfa_diff_birth_inf12 		-0.19 -0.45 – 0.07 0.148 0.09 -0.85 – 1.04 0.846 -0.51 -1.05 – 0.04 0.068 -1.18 -2.84 – 0.48 0.161 -0.41 -0.98 – 0.15 0.153 -1.04 -2.70 – 0.62 0.217 0.06 -0.28 – 0.40 0.747 0.25 -0.97 – 1.47 0.687
Observations 370 100 370 100 370 100 370 100
R2 / R2 adjusted 0.196 / 0.188 0.021 / -0.010 0.088 / 0.078 0.032 / 0.001 0.023 / 0.012 0.023 / -0.007 0.004 / -0.007 0.006 / -0.025 
#weight 83-91 interaction models 
par(mfrow=c(2,2))
plot(grim_weight_83_91_intxn_f)

par(mfrow=c(2,2))
plot(grim_weight_83_91_intxn_m)

par(mfrow=c(2,2))
plot(pheno_weight_83_91_intxn_f)

par(mfrow=c(2,2))
plot(pheno_weight_83_91_intxn_m)

par(mfrow=c(2,2))
plot(han_weight_83_91_intxn_f)

par(mfrow=c(2,2))
plot(han_weight_83_91_intxn_m)

Looked at another way (weight and height, females only, 83-91)

sjPlot::tab_model(grim_weight_83_91_intxn_f, 
                  pheno_weight_83_91_intxn_f, 
                  han_weight_83_91_intxn_f,
                  horv_weight_83_91_intxn_f,
                  grim_height_83_91_intxn_f, 
                  pheno_height_83_91_intxn_f, 
                  han_height_83_91_intxn_f, 
                  horv_height_83_91_intxn_f)
  AgeAccelGrim AgeAccelPheno EEAA IEAA AgeAccelGrim AgeAccelPheno EEAA IEAA
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) -1.07 -1.85 – -0.28 0.008 -0.26 -1.89 – 1.38 0.757 0.21 -1.49 – 1.91 0.807 0.50 -0.52 – 1.53 0.335 -0.92 -1.41 – -0.43 <0.001 -0.18 -1.20 – 0.84 0.730 -0.11 -1.17 – 0.96 0.844 -0.04 -0.67 – 0.60 0.904
wfa_diff_inf12_91 0.66 -0.04 – 1.35 0.064 1.47 0.03 – 2.92 0.046 1.39 -0.12 – 2.89 0.072 -0.26 -1.17 – 0.64 0.568
wfa_diff_birth_inf12 0.07 -0.23 – 0.36 0.650 -0.05 -0.66 – 0.56 0.872 -0.27 -0.91 – 0.37 0.414 -0.23 -0.61 – 0.16 0.245
was_preg_no_na [Yes] 2.79 2.19 – 3.38 <0.001 3.49 2.25 – 4.72 <0.001 1.00 -0.29 – 2.28 0.129 0.02 -0.76 – 0.80 0.958 2.76 2.17 – 3.35 <0.001 3.42 2.18 – 4.66 <0.001 0.86 -0.43 – 2.15 0.192 0.01 -0.76 – 0.78 0.974
wfa_diff_inf12_91 *
wfa_diff_birth_inf12 		-0.19 -0.45 – 0.07 0.148 -0.51 -1.05 – 0.04 0.068 -0.41 -0.98 – 0.15 0.153 0.06 -0.28 – 0.40 0.747
hfa_diff_inf12_91 0.25 -0.27 – 0.76 0.344 0.02 -1.06 – 1.09 0.975 0.24 -0.88 – 1.36 0.671 -0.03 -0.70 – 0.63 0.919
hfa_diff_birth_inf12 -0.01 -0.25 – 0.23 0.926 0.08 -0.42 – 0.58 0.756 0.25 -0.26 – 0.77 0.335 -0.07 -0.38 – 0.24 0.673
hfa_diff_inf12_91 *
hfa_diff_birth_inf12 		0.04 -0.16 – 0.24 0.680 -0.07 -0.49 – 0.34 0.723 -0.13 -0.55 – 0.30 0.566 0.12 -0.13 – 0.38 0.339
Observations 370 370 370 370 370 370 370 370
R2 / R2 adjusted 0.196 / 0.188 0.088 / 0.078 0.023 / 0.012 0.004 / -0.007 0.192 / 0.183 0.077 / 0.067 0.014 / 0.003 0.012 / 0.001 
vif(grim_height_83_91_f)
##    hfa_diff_inf12_91 hfa_diff_birth_inf12       was_preg_no_na 
##             1.183167             1.179847             1.003111
vif(grim_height_83_91_m)
##    hfa_diff_inf12_91 hfa_diff_birth_inf12 
##             1.203573             1.203573
vif(pheno_height_83_91_f)
##    hfa_diff_inf12_91 hfa_diff_birth_inf12       was_preg_no_na 
##             1.183167             1.179847             1.003111
vif(pheno_height_83_91_m)
##    hfa_diff_inf12_91 hfa_diff_birth_inf12 
##             1.203573             1.203573
vif(han_height_83_91_f)
##    hfa_diff_inf12_91 hfa_diff_birth_inf12       was_preg_no_na 
##             1.183167             1.179847             1.003111
vif(han_height_83_91_m)
##    hfa_diff_inf12_91 hfa_diff_birth_inf12 
##             1.203573             1.203573
vif(grim_height_91_02_f)
##       hfa_diff_91_02 hfa_diff_birth_inf12       was_preg_no_na 
##             1.056932             1.051030             1.006030
vif(grim_height_91_02_m)
##       hfa_diff_91_02 hfa_diff_birth_inf12 
##             1.020674             1.020674
vif(pheno_height_91_02_f)
##       hfa_diff_91_02 hfa_diff_birth_inf12       was_preg_no_na 
##             1.056932             1.051030             1.006030
vif(pheno_height_91_02_m)
##       hfa_diff_91_02 hfa_diff_birth_inf12 
##             1.020674             1.020674
vif(han_height_91_02_f)
##       hfa_diff_91_02 hfa_diff_birth_inf12       was_preg_no_na 
##             1.056932             1.051030             1.006030
vif(han_height_91_02_m)
##       hfa_diff_91_02 hfa_diff_birth_inf12 
##             1.020674             1.020674
vif(grim_weight_83_91_f)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12       was_preg_no_na 
##             1.275654             1.284605             1.008342
vif(grim_weight_83_91_m)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12 
##             1.053904             1.053904
vif(pheno_weight_83_91_f)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12       was_preg_no_na 
##             1.275654             1.284605             1.008342
vif(pheno_weight_83_91_m)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12 
##             1.053904             1.053904
vif(han_weight_83_91_f)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12       was_preg_no_na 
##             1.275654             1.284605             1.008342
vif(han_weight_83_91_m)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12 
##             1.053904             1.053904
vif(horv_weight_83_91_f)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12       was_preg_no_na 
##             1.275654             1.284605             1.008342
vif(horv_weight_83_91_m)
##    wfa_diff_inf12_91 wfa_diff_birth_inf12 
##             1.053904             1.053904
vif(grim_height_83_91_intxn_f)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               4.055098                               1.561438 
##                         was_preg_no_na hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               1.003588                               5.066738
vif(grim_height_83_91_intxn_m)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               3.141898                               1.471774 
## hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               3.799584
vif(pheno_height_83_91_intxn_f)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               4.055098                               1.561438 
##                         was_preg_no_na hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               1.003588                               5.066738
vif(pheno_height_83_91_intxn_m)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               3.141898                               1.471774 
## hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               3.799584
vif(han_height_83_91_intxn_f)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               4.055098                               1.561438 
##                         was_preg_no_na hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               1.003588                               5.066738
vif(han_height_83_91_intxn_m)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               3.141898                               1.471774 
## hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               3.799584
vif(horv_height_83_91_intxn_f)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               4.055098                               1.561438 
##                         was_preg_no_na hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               1.003588                               5.066738
vif(horv_height_83_91_intxn_m)
##                      hfa_diff_inf12_91                   hfa_diff_birth_inf12 
##                               3.141898                               1.471774 
## hfa_diff_inf12_91:hfa_diff_birth_inf12 
##                               3.799584
vif(grim_height_91_02_intxn_f)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.111082                            1.176695 
##                      was_preg_no_na hfa_diff_91_02:hfa_diff_birth_inf12 
##                            1.008046                            4.442128
vif(grim_height_91_02_intxn_m)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.117563                            1.257104 
## hfa_diff_91_02:hfa_diff_birth_inf12 
##                            4.576882
vif(pheno_height_91_02_intxn_f)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.111082                            1.176695 
##                      was_preg_no_na hfa_diff_91_02:hfa_diff_birth_inf12 
##                            1.008046                            4.442128
vif(pheno_height_91_02_intxn_m)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.117563                            1.257104 
## hfa_diff_91_02:hfa_diff_birth_inf12 
##                            4.576882
vif(han_height_91_02_intxn_f)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.111082                            1.176695 
##                      was_preg_no_na hfa_diff_91_02:hfa_diff_birth_inf12 
##                            1.008046                            4.442128
vif(han_height_91_02_intxn_m)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.117563                            1.257104 
## hfa_diff_91_02:hfa_diff_birth_inf12 
##                            4.576882
vif(horv_height_91_02_intxn_f)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.111082                            1.176695 
##                      was_preg_no_na hfa_diff_91_02:hfa_diff_birth_inf12 
##                            1.008046                            4.442128
vif(horv_height_91_02_intxn_m)
##                      hfa_diff_91_02                hfa_diff_birth_inf12 
##                            4.117563                            1.257104 
## hfa_diff_91_02:hfa_diff_birth_inf12 
##                            4.576882
vif(grim_weight_83_91_intxn_f)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                               5.939198                               1.311235 
##                         was_preg_no_na wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                               1.008646                               5.363499
vif(grim_weight_83_91_intxn_m)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                              12.384711                               1.233082 
## wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                              13.154471
vif(pheno_weight_83_91_intxn_f)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                               5.939198                               1.311235 
##                         was_preg_no_na wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                               1.008646                               5.363499
vif(pheno_weight_83_91_intxn_m)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                              12.384711                               1.233082 
## wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                              13.154471
vif(han_weight_83_91_intxn_f)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                               5.939198                               1.311235 
##                         was_preg_no_na wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                               1.008646                               5.363499
vif(han_weight_83_91_intxn_m)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                              12.384711                               1.233082 
## wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                              13.154471
vif(horv_weight_83_91_intxn_f)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                               5.939198                               1.311235 
##                         was_preg_no_na wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                               1.008646                               5.363499
vif(horv_weight_83_91_intxn_m)
##                      wfa_diff_inf12_91                   wfa_diff_birth_inf12 
##                              12.384711                               1.233082 
## wfa_diff_inf12_91:wfa_diff_birth_inf12 
##                              13.154471